The Detectability of Exo-Earths and Super-Earths Via Resonant 

Signatures in Exozodiacal Clouds 



Christopher C. Stark^ and Marc J. Kuchner^ 
ABSTRACT 

Directly imaging extrasolar terrestrial planets necessarily means contending 
with the astrophysical noise of exozodiacal dust and the resonant structures cre- 
ated by these planets in exozodiacal clouds. Using a custom tailored hybrid 
symplectic integrator we have constructed 120 models of resonant structures cre- 
ated by exo-Earths and super-Earths on circular orbits interacting with colli- 
sionless steady-state dust clouds around a Sun-like star. Our models include 
enough particles to overcome the limitations of previous simulations that were 
often dominated by a handful of long-lived particles, allowing us to quantitatively 
study the contrast of the resulting ring structures. We found that in the case 
of a planet on a circular orbit, for a given star and dust source distribution, the 
morphology and contrast of the resonant structures depend on only two param- 
eters: planet mass and ^Ja^j(5^ where Op is the planet's semi-major axis and /3 
is the ratio of radiation pressure force to gravitational force on a grain. We con- 
structed multiple-grain-size models of 25,000 particles each and showed that in a 
collisionless cloud, a Dohnanyi crushing law yields a resonant ring whose optical 
depth is dominated by the largest grains in the distribution, not the smallest. 
We used these models to estimate the mass of the lowest-mass planet that can 
be detected through observations of a resonant ring for a variety of assumptions 
about the dust cloud and the planet's orbit. Our simulations suggest that planets 
with mass as small as a few times Mar's mass may produce detectable signatures 
in debris disks at Op > 10 AU. 

Subject headings: catalogs — circumstellar matter — infrared: stars — inter- 
planetary medium — methods: N-body simulations — planetary systems 



^Department of Physics, University of Maryland, Box 197, 082 Regents Drive, College Park, MD 20742- 
4111, USA; starkc@umd.edu 

^NASA Goddard Space Flight Center, Exoplanets and Stellar Astrophysics Laboratory, Code 667, Green- 
belt, MD 20771 



- 2 - 



Introduction 



A number of proposed experiments like the Terrestrial Planet Find er (TPF) aim to 



direc tly image the scattered and emitted hght from extrasolar planets (iLawson fc Traub 



Agol 


2007; 


Beckwith 


2007) 



analogous to zodiacal dust in our solar system (e.g. 
dust in the solar system consists of ~ 1 — 10 fim dust grains rele ased through asteroidal 



collisions and the outgassing of comets (e.g. ISchramm et al.lll989l). T his dust forms the 



zodi acal cloud, extendin g from the solar corona (e.g. 
(e.g. iKriiger et al.l 119991 ). 



Mann et al 



20001 ) to beyond Jupiter 



Our zodiacal cloud exhibits several structures in terpreted as dynamical signatures of 
planets (jPermott et al.l Il985l . 1 1994 : iReach et al.l Il995l). Several dusty d isks around nearby 



main-sequence st ars show similar structures (e.g. iGreaves et al.l Il998l : IWilner et al.l 12002 
Kalas et al.ll2005l ). This trend suggests that exozodiacal clouds may be full of rings, clumps 
and other asymmetries caused by planets and other phenomena. 

This situation raises some important questions. Will the structures in e xozodiacal clouds 



be harmful astrophys ical noise for direct imaging of extrasolar planets (IBeichmanl Il996 



Beichman et al.lll999l )? Or can the dyna mical signatures of planet s in these clouds help us 
find otherwise undetectable planets (e.g. iKuchner fc Holmanll2003l )? 



Sever al studies have exarained the geometry of res onant signatures of planets in debris 



disks (e.g. iKuchner &: Holmanll2003l : iReche et al.ll2008l ). However, most simulations cannot 



quantitatively study the contrast in these structures: how bright they are relative to the 
background cloud. We need to model the contrast of the structures in exozodiacal clouds to 
understand their roles as astrophysical noise and as signposts of hidden planets. However, 
accurately simulating the contrast of thes e structures demands corn putational resources that 
have only recently become available (e.g. iDeller fc Maddisonll2005l ). 



In this paper we examine the contrast of resonant structures induced by planets in 
steady-state exozodiacal clouds and the detectabihty of these structures via direct imaging. 
We simulate high-fidelity images of coUisionless exozodiacal clouds containing a terrestrial- 
mass planet — an exo-Earth or super-Earth. By using roughly an order of magnitude more 
particles than most previous simulations, we overcome the Poisson noise associated with 
constructing histograms of the column density and populating the external mean motion 
resonances (MMRs) of planets. We use our simulations to estimate the minimum planet 
mass that can be indirectly detected via observations of these structures as a function of 
the planet semi-major axis and dominant grain size under the assumption of circular planet 
orbits. Our models apply to exozodiacal clouds less than a few hundred times the optical 
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depth of the solar zodiacal cloud, clouds for which the collision time is shorter than the 
Poynting- Robert son (PR) time for typical grains. 

Section 2 of this paper describes our numerical techniques. We present a synthetic 
catalog of resonant debris disk structures in Section 3. We describe our multiple-particle- 
size cloud models and discuss their detectability in Section 4. In Section 5, we discuss the 
limitations of our simulations; we summarize our conclusions in Section 6. 



2. Numerical Method 



Dust grains in the inner solar system are primarily released from parent bodies via colli- 
sions or outgassing. Radiation pressure ejects the smallest particles from the sola r system in 
a dynamical time while the larger particles slowly spiral inward due to PR drag (IRobertson 
19371 : iBurns et al.lll979l ). During their spiral toward the Sun, particles may become tem- 
pora rily trapped in the MM Rs of planets, extending their lifetimes by a factor of a few to 
ten (j Jackson &: ZooM Il989l ). This trapping locally enhances the particle density, creating 
structures within the zodiacal cloud , which have been described as circumsolar rings, bands, 
and clumps (e.g. iKelsall et al.lll998l ). 



To model these types of structures in exozodiacal clouds we numerically integrated the 
equation of motion of dust particles. The equation of motion f or a perfectly ab sorbing 
particle orbiting a star of mass is given to first order in w/c by iRobertsonI (119371 ): 



sw)/3 Gm^ 



rr 



:i^ 



where r and v are the heliocentric position and velocity of the parti cle and sw is the ratio 
of solar wind drag to PR drag. We assume a value for sw of 0.35 (iGustafsonI Il994j ) . For 
perfectly absorbing spherical particles in the vicinity of the Sun, /3 ~ 0.57 /ps, where p is the 
mass density of the particle in g cm~^ and s is the radius in pm. 



2.1. A Customized Hybrid Symplectic Integrator 



We imple mented a c ustorn ized hybrid symplectic i ntegr ator to perform our numerical 
integrations. I Chambers! (119991 ) . hereby referred to as IC99l . introduced hybrid symplectic 
integration as a method for dealing with close encounters in an efficient n-body code. Sym- 
plectic integrators rely on splitting the Hamiltonian into two easily integrable portions — a 
dominant term, Hd, and a smaller perturbative term. Hp. However, in the n-body problem. 
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Hp may exceed i^D during close encounters. Hybrid symplectic integrators overcome this 
problem by effectively switching from a symplectic integrator to an alternate integrator (e.g. 
Buhrsch-Stoer). 

The hybrid method reduces the perturbative term of the Hamiltonian, Hp, by a factor 
Kijij), where r^j is the distance between the two bodies in question, to ensure that the 
perturbative term remains relatively small. The integrator includes the remaining portion of 
the perturbative term, ^. Hpij[l — ii"(rjj)], in the dominant term which is then integrated 



using a method of choice. The "changeover function, 
varies from for rtj < r^it to unity for r^j > rcrit- 



K(rij), is a smooth function that 



Using a hybrid integrator requires c hoosin g a changeover function and a value for rcrit- 
We use the same changeover function as IC99l . We assign a different value of rcrit,i to each 
body, calculated as the larger of 3i?H,i and rwj, where i?H,j and Vi are the Hill radius and 
velocity of the i^^ body, respectively, and r is the time step of the integrator. We then 
calculate the critical distance for a pair of bodies as rcnt,ij = rcrit,i + '^critj- 

Our integrator also incorporates the effects of radiation pressure, PR drag, and solar 
wind drag. We implement radiation pressure as a correction to the effective stellar mass 
(c f. Eg. 1) and treat the drag effe cts as an additional te rm in Hp , in much the same way 



as 



Moro-Martm fc Malhotral ((2002), hereby referred to as IMMM02I . We also use democratic 



heliocentric (DH) coordinates, composed of the barycentric momenta and heliocentric posi- 
tions, because of their relative ease of implementation. This choice introduces an additional 
perturbativ e term to the Ham iltonian due to the motion of the star with respect to the 
bary center (IDuncan et al.lll998l ). 



2.2. Comparison of Integrator with Previous Results 

We checked our integrator using a variety of standard tests. We checked the energy and 
Jacobi constant conservation with the drag terms turned off and examined the evolution of 
dust particles' orbital elements under our implementation of drag effects. We also compared 
our hybrid integrator to a Bulirsch-Stoer integrator by examining the path of an individual 
test particle during a close encounter and by examining the statistics of a cloud of particles 
in a collisionless disk containing a planet. 

We tested energy conservation in our integration code by integrating the orbits of the 
four outer planets and the Sun for 3 x 10^ years using a time step of 0.15 years. The energy 
error was bounded with a mean value of AE/E ^ 3 x 10~^. 
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Duncan et al.l (Il998l ). which we will refer to as IDLL98I . tested the relative conservation 
of energy in their symplectic integrator as a function of planet perihelion distance. We 
replicated their tests using our code. Figured] shows the relative energy error in an integration 
of the orbit of Jupiter for 3 x 10^ years using a time step of 0.15 years. We initially placed 
Jupiter at aphehon. Figure [1] also shows the results of integrating the orbits of Jupiter 
and Saturn under the same conditions. With the DH method, the perturbative solar term 
increases as the perihelion distance of the planet decreases, causing the fractional energy 
error to incr ease sim ilarly. The fractional energy errors shown in Figure [T] agree with those 



obtained by IDLL98 



We checked the conservation of the Jacobi constant by integratin g particles in the Sun- 
Neptune system. We found results consistent with those of IMMM02I . Particles that did not 
undergo close encounters conserved the Jacobi constant at the level of ~ 10^^ to 10^^. 



To test our implementation of PR drag, we replicated a test performed by IMMM02 . 
We integrated the orbit of a particle with (3 = 0.2 and sw = 0.35 in the presence of the 
Sun. Figure E] show s the semi -major axis and eccentricity as functi ons of time. These result s 
match the results of lMMM02l and agree with the analytic solution (IWyatt fc Whipplelll950l ). 



We tested the performance of our hybrid scheme by integrating the orbit of comet 

P / Ote rma in a c lose encounter with Jupiter, which has been done previously bv lMichel fc Valsecchil 
( 1996 ) and|C99|. The initial conditions for both bodies can be found in Table 3 of lMichel fc Valsecchi 
( 119961 ) . Figure [3] shows the path of comet P / Oterma for several values of integration time 
step r as seen in the frame co-rotating with Jupiter, which is located at the origin. These 
results are similar to those obtained by IC99l . Our code shows a minor improvement over the 
other codes, most noticeable in the r = 100 d ays ca se, that is likely only due to differences 
in the calculation of the changeover distance. IC99l explicitly sets rcrit = 3i?H for this test; 
we used our prescription for rcrit as described in Section 12.11 



2.3. Test Simulations of a Steady-State Exozodiacal Cloud 



We directly compared simulations of resonant structures made with our hybrid integra- 
tor to simulations made with a Bulirsch-Stoer integrator. During the integrations we recorded 
the coordinates of each particle in a 2-D histogram at regular intervals. This histogram mod- 
els the surface density distribution in a steady-state cloud. Since we only modeled planets 
on circular orbits, we simply recorded the coordinates in the frame co-rotating with the 



planet. This techii i que has been widely used by dust cloud modelers ( Dermott et al 



Liou fc Zookl Il999l : iMoro-Martm fc Malhotral l2002l : IWilner et all 12002 



11994 



Deller &: Maddison 
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20051). 

Figure Hlshows two histograms, one for each integrator, for simulations of 1,000 particles 
each in the presence of the Sun/Earth system. We used a histogram bin width of 0.0175 AU. 
For these simulations we chose (5 = 0.02 and initially released the particles with semi-major 
axis, Odust, distributed uniformly between 3 and 5 AU, eccentricity, e, uniformly distributed 
from 0.0 to 0.1, and inclination, i, uniformly distributed between 0° and 6°. We used a 
symplectic time step of 0.02 years and recorded the particle locations every 250 years. 

Except for a small number of pixels, the middle panel of the figure (simulation using the 
hybrid integrator) looks qualitatively very similar to the left-most panel (simulation using 
the Bulirsch-Stoer integrator). The right panel of Figure H] shows the difference of these two 
images divided by the y/n Poisson noise expected for each pixel where n is the number of 
particles in the pixel. This figure demonstrates that the differences between the two models 
are nearly consistent with the Poisson noise of the histograms. The two integrators resulted 
in histograms with minor structural differences, but the hybrid symplectic integrator runs a 
few times faster. 

Besides pixel-to-pixel Poisson noise in th e histogram, this method is also sensitive to 



noise in the population of MMRs. IMMM02I showed that the population of the dominant 



resonances varied by a factor of ~ 3 among sets of 100-particle simulations of Kuiper Belt 
dust interacting with Neptune. This noise probably causes the differences between the two 
simulations shown in Figure H] beyond those attributable to pixel-to-pixel Poisson noise. 
Although simulations of 100 particles may acquaint us with the generic geometry of debris 
disk structures, we cannot use them to predict ring contrasts; to model the contrast in a 
resonant cloud feature we must include enough particles to accurately populate the MMRs. 

We solved this problem by using more particles. We used the 420-processor Thunderhead 
cluster at NASA Goddard Space Flight Center to perform simulations of 5,000 particles each. 
Figure [5] shows the population of MMRs for three independent 5,000- particle si mulations of 
the Sun and four outer planets using the same initial conditions as IMMM02I . Simulating 
5,000 particles reduced the difference between MMR populations for the three simulations 
to less than 7% for the dominant 2:1 and 3:2 MMRs, allowing us to synthesize high-fidelity 
images and quantitatively study the resonant ring structures. 
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Simulations & Results 



3.1. Cataloging Debris Disk Structure 



To explore the range of different types of structures formed by terrestrial-mass planets, 
we performed 120 simulations of dust interacting with single planets on circular orbits. The 
simulations used 5,000 particles each and covered six values of planet mass, Mp (0.1, 0.25, 
0.5, 1.0, 2.0, and 5.0 M®), four values of planet semi-major axis, Op (1, 3, 6, and 10 AU), and 
five values of (3 (0.0023, 0.0073, 0.023, 0.073, and 0.23) corresponding to spherical sihcate 
particles ranging in radius from ~ 1 — 120/im. We released the particles on orbits with 
semi-major axes uniformly distributed between 3.5 and 4.5 times the semi-major axis of the 
planet's orbit — well outside of the strongest MMRs. We used initial eccentricities uniformly 
distributed between and 0.2, initial inclinations uniformly distributed between and 20°, 
and the longitude of the ascending node, Q, and the argument of pericenter, u, uniformly 
distributed between and 27r. We considered planet semi-major axes of 1 to 10 AU because 
typical designs for TPF can detect an exozodiacal cloud with 10 times the optical depth of 
the solar zodiacal cloud over roughly that range of circumstellar radii (ILevine et al.ll2006l ). 



We chose these initial conditions to model only dynamically-cold dust, i.e. Cdust ^ 0.2 
and idust ^ 20°, since this component of a dust cloud is the dominant contributor to resonant 
ring structure. We neglect dynam ically hot dust wit h the idea that it can always be added 



in later as a smooth background (iMoran et al.ll2004l ). The asteroid belt probably produces 



Liou et al. 


1995; 


Ipatov et al. 


2007) 



more dynamically hot cloud component (e.g. 
only steady-state dust clouds, assuming dust is continually replenished, and ignore transient 
coUisional events. 

Figure M shows some examples of the histograms from our simulations, which reveal a 
wide range of trapping behavior. Some histograms show no azimuthal or radial structure, 
wh ile others show high contra st rings. All of the patterns are Type I structures as identified 
by lKuchner fc Holmaru (120031 ). 



Several general trends emerged. The ring contrast increased with increasing planet 
mass. Reducing (3 also enhanced trapping, as did increasing the planet's semi-major axis. 

These last two trends can be explained by comparing the libration time of a given MMR 
to the PR time. The PR time scales as ctdust//^' while the libration time for a given resonance 
scales as a^usf where Odust is the semi-major axis of a dust grain's orbit. The ratio of these 
quantities yields y/adust/P-, a pa rameter that measures the degree to which resonant trapping 
is adiabatic (e.g. lHenrardlll982l ): the trapping becomes more adiabatic and more efficient at 
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greater distances from the star, and for larger particles. We discuss this phenomenon further 
in Section [331 below. 

In addition to following these trends, all of the simulated ring structures, like those 
shown in Figure [6l share some salient features: 

1. For cases in which even a modest amount of trapping occurs (azimuthally averaged 
contrasts ^ 1.3 : 1), the ring structures exhibit a sharp inner edge at ~ 0.83ap. This 
feature probably appears because the eccentricities of particles trapped in exterior 
MMRs are typically pumped up to a limiting value before a close encounter with the 
planet ejects them from resonance. For a particular MMR, all particles, regardless of 
f3, tend to ap proach the same limitirig ecce ntricity and accordingly, a similar pericen- 
ter distance (iBeauge &: Ferraz-Melldll994l ). The limiting eccentricities are such that 
the limiting pericenter distances are nearly equal for the dominant resonances (e.g. 
the 2:1 and 3:2 resonances have limiting pericenter distances of 0.823ap and 0.827ap, 
respectively), creating the ring structure's sharp inner edge. 

2. A gap in the ring structure, a local minimum in the surface density, appears around 
the planet. If we define gap width as the FWHM of the minimum in the azimuthal 



surface density profile at r 



we find that the gap width is linearly proportional 



to the contrast of the ring, as shown in the left panel of Figure [71 A linear fit to the 
data shown in this figure gives Wgap ~ 10° x Caa,ie for Caa.ie > 1-6, where Wg^p is the 
gap width in degrees and Caa,ie is the azimuthally averaged inner-edge contrast (see 
Section [32D- 

3. The rings show a leading-trailing asymmetry. The trailing side of the ring structure is 
noticeably denser than the leading side, and the structure is rotationall y shifted in the 



progr ade direction causing the trailing side to be closer to the planet (jPermott et al. 



19941 ). To examine the leading-trailing asymmetry caused by a prograde shift of the 



ring structure, we measured the azimuthal offset of the center of the gap described 
above from the planet. The right panel of Figure [7] shows these measured prograde 
shifts of our simulations. iKuchner fc HolmanI (120031 ) showed for a particular first order 
exterior MMR, 



sm (po oc 



(2) 



where 0o is the prograde shift of the pericenter. Therefore, we plotted the sine of each 
measured prograde shift against /3(1 — (3) / {Mp^Ja^) in the right panel of Figure [3 Our 
data reveal the approximate proportionality 

> 1 0.5 



sm I 



'rmg 



OC 



(3) 
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where (f)j:mg is the measured prograde shift of the ring structure. While the relationship 
in Equation [2] holds for a single MMR, it does not strictly apply to a given ring 
structure which consists of several well-populated MMRs. The relative populations of 
these MMRs are also functions of Mp, a^, and f3. However, this situation seems to 
preserve a power-law relationship between sin^nng and /3 (1 — /?) / (Mpy/a^), as shown 
in Equation [31 

4. The radial width of the ring increases with the contrast of the ring, ranging from a few 
percent of Op to ~ l.Gop in the highest contrast case. As the trapping probabilities of 
all the MMRs increase, MMRs farther from the planet's orbit become populated. For 
this reason, the outer-edge of the ring structure differs significantly among simulations. 
The outer-edge can be quite blurry or very well-defined, making the radial width of a 
ring structure difficult to quantify. 



Our catalog of debris disk structures induced by terrestrial-mass planets is publicly avail- 



able online at http: / /asd. gsfc.nasa.gov/Christopher.Stark/ catalog. php| This online catalog 
also contains images synthesized from the density distributions in scattered light and 10 /im 
thermal emission assuming blackbody grains. Future studies of resonant ring structures with 
TPF or other experiments can use our catalog to interpret dust cloud patterns in terms of 
planet and dust parameters, assuming the observed image is dominated by a single grain 
size. We envision the followin g process, inspired by recent papers on disks observed with the 
Hubble Space Telescope (e.g. IClampin et al.ll2003l : iKalas et al.ll2005l ): 



1. Deproject the image to remove inclination effects. 

2. Remove any smooth backgrounds by a power law fit. 

3. Estimate the dominant grain size in the resonant ring using infrared photometry or 
other methods. 



4. Compare the image of the disk to the online catalog to constrain the planet's mass and 
location. 



3.2. Ring Contrast 

We considered three different metrics for describing the ring contrast in our simulations: 

Civiax: The surface density of the ring at its densest point divided by the surface density 
of the background cloud 



10 



CAA,Max: The maximum value of the azimuthally averaged surface density divided by 
the surface density of the background cloud 

Caa,ie: The azimuthally averaged surface density at the inner edge of the ring divided 
by the surface density of the background cloud 

We calculated the above contrast metrics for all 120 simulations. We measured the 
surface density of the background cloud at a circumstellar distance r ^ 0.8ap. The surface 
density of the background cloud was nearly constant inside and outside of the ring, but did 
exhibit a small local minimum near r ^ O.Sop in a few cases. 

To calculate Cmsx, we must search for the densest pixel, which introduces a bias toward 
pixels that exhibit an extreme amount of Poisson noise. To reduce this noise we averaged the 
surface density over nine pixels centered on the densest point. Using CAA.Max or Caa.ie, on 
the other hand, automatically averages over the effects of Poisson noise in our simulations. 

Figure [8] shows two examples of how the contrast, Caa,ie, depends on planet mass and /?. 
Both plots show a similar behavior with three distinct regions: a no-trapping regime (contrast 
~1), a transitional regime, and a saturation regime (maximum contrast). The saturation 
regime is of particular significance. Our results suggest that within the range of parameters 
investigated, for a given value of (3, all contrasts converge to the same value for large planet 
masses independent of planet semi-major axis, i.e. the contrast becomes "saturated" and 
increasing the planet's semi-major axis has little effect on the contrast. The right panel in 
Figure [8] illustrates this behavior; all four contrast curves, each of which corresponds to a 
different planet semi-major axis, approach the same value of ~ 7 near Mp = 5 M^. Similarly, 
for a given planet mass, contrasts converge to the same value for small [3 independent of Op, 
as shown in the left panel in Figure [HI The morphology of the structure can vary, but the 
contrast of the ring structure is roughly constant in these saturation regimes. 



As we mentioned above, dividing the PR time by the libration time of a given MMR 
yields a parameter, ^ya^/j3, that indicates the degree to which the resonant trapping is 




adiabatic. We plot the contrast in our simulations as a function of this parameter in Figure 
O This figure demonstrates that for Mp < 5 and for a given distribution of parent body 
orbital elements, the ring contrast is a function of only two parameters: planet mass and 



The morphology of the resonant rings is also, to good approximation, a function of only 



3.3. Adiabaticity 
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the planet mass and ^Ja^/ (3. The models shown in the two right panels of Figure [6] illustrate 
this phenomenon; for both models yfa^/(3 ~ 137 AU^^^. These two models have the same 
morphology to a level consistent with pixel-to-pixel Poisson noise. Note that for small (3 in 
Equation [31 the prograde shift is approximately a function of (y^//5) ^ for a given planet 
mass. 

For large values of and Mp, the morphology and contrast of the ring structures are not 
simple functions of ^Ja^jfi. Simulations with large values of /3, but equal values of ^/a^/[3 
(e.g. Vl AU/0.073 and VlO AU/0.23) show morphological differences, including differences 
in prograde shift. Our simulations with Mp = 5 also show contrast differences among 
rings with equal values of Ja^jfi. 



WyattI (I2OO3I ) investigated resonant trapping in MMRs for a syste r n of p lanetesimals 



exterior to an outward migrating planet on a circular orbit. IWyattI (120031 ) plotted the 



trapping probability for a single MMR in his model as a function of migration rate and planet 
mass and found it could be well approximated by a function of the form P = [1 + (dp/pi)^^]"^, 
where dp is the migration rate and the parameters pi and p2 are power laws in planet mass. 
Our trapping scenario assumes dust migrating inward toward the planet, but the concept is 
similar. Since contrast is closely related to trapping probability, we decided to fit the data 
shown in Figure [9] with a function of the form 



C = 1 + pi 1 + 



P2 



Pa 



-1 



(4) 



inspired by IWyattI (120031 ). Each of the three parameters, Pi, is a power law in planet mass 
of the form pi = pi^iMp''^ . We fit all 120 contrast measurements with this six-parameter 
function for each of the three contrast metrics. The best fits, two of which are shown in 
Figure [9], are: 



Caa,ie: Pi ^ 

CAA.Max^ Pi 



x0.19 



4.38 (Mp/Me)-",P2 



207 {Mp/M^ 



^1.17 



4.54 (Mp/Me)°-'', P2 ^ 205 (Mp/Mg 



-1.17 



P3 ^ 2.05 (Mp/Mgj) 
P3~1.63 (Mp/Mg 



0.11 



n0.19 



CMax: Vi ~ 6.23 (Mp/Me)°-'', P2 ^ 164 {MJM, 



-1.09 



P3~1.72 (Mp/Mg 



n0.05 



Equation [D, combined with the above values, summarizes our results for all combinations 
of planet mass, planet semi-major axis and (3 we simulated. Figure [9] shows the inner- 
edge contrast, Caa.ie, deviates significantly from the fits for large Mp and large yfa^/jS. 
The increased trapping efficiency for MMRs with these massive planets likely enhances the 
population of MMRs farther from the planet's orbit and depletes the inner MMRs that cause 
the sharp inner edge. 
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Multi-Particle-Size Models 



4.1. Composite Simulations 

We used our 120 simulations to produce 20 multiple-particle-size dust cloud models 
by fo rming weighted sums of the histograms assuming a Dohnanyi distribution of particle 
sizes (jPohnanyil Il969l ) . Each of th ese compos i te mo dels effectively utilizes 25,000 particles. 
Exactly how we apply the ideas in iDohnanyil (119691 ) has profound effects on our composite 
models, so we present here two different kinds of models. 

First, we assembled a composite model in which the particles are initially released from 
their parent bodies according to a crushing law, and do not undergo any further collisional 
processing as they spiral inward. This scenario models a sparse disk with a belt of dust- 
producing material, like our own zodiacal cloud. The crushing law for ast e roid m aterial at 
micron sizes is unknown, so we choose the crushing law used by IDohnanyil (119691 ): 



dN 
ds 



oc s 



(5) 



where dN is the number of particles with radius s in a bin of width ds, and a = 3.4. 



We calculate the optical depth, r, for our composite models from r = J2i ua^i^i, where 
Wi, Ai, and (jj are the weighting factor, particle cross-section, and surface number density 
of the single-particle simulation, respectively. We assume that the cross-section of each 
particle is Ai oc The crushing law in Equation |5] implies a weighting factor for the i^^ 
histogram of Wi = /?f ~^A/3j, where A/3j is the width of the bin. For a constant logarithmic 
spacing in j3, like the spacing we used in our simulations, and the Dohnanyi crushing law. 

Larger particles have longer PR times, so in the absence of collisional processing, their 
density is enhanced by a factor of under our assumption of a steady-state cloud model. 
One might expect that this effect must be included in the weighting factor. However, our 
simulations include this effect automatically as long as we keep the frequency with which 
particle locations are recorded constant among all of our simulations. We did, in fact, vary 
the recording frequency with the PR time, but we corrected for the differences in recording 
frequency before summing the histograms. 

Figure Uni shows the optical depth of one of our 20 composite models (Mp = 2.0 M^, 
Op = 6.0 AU), together with the optical depths of single-particle-size models using only the 
smallest and largest par t icle s izes included in the composite model. Although the crushing 
law used by IDohnanyil (Il969l) favors small er particles by number, even more than some 
empirical crushing laws (IDurda et al.ll2007l ). the optical depth in the composite models is 
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dominated by the largest particles. This situation occurs because the larger particles are 
both longer lived (tpR oc and more likely to be trapped in MMRs. Hence, the upper 
left panel of the figure closely resembles the lower right panel. 

Next, for the purpose of illustration, we ignored the initial size distribution of dust 
particles and forced the disk to obey a size distribution of 

dN 

— oc s-^-^ 6 
as 

at a radius of ~ 3ap from the star. This scenario probably doesn't have a physical inter- 
pretation, but it illustrates an interesting phenomenon: how resonant trapping tends to sort 
particles by size. We enforce the size distribution at one location within the disk, but the 
size distribution will not follow a Dohnanyi distribution elsewhere in the disk. 

The top right panel in Figure [10] shows the optical depth of an example of this kind of 
composite cloud, normalized to a Dohnanyi distribution at ~ 18 AU. Models constructed in 
this fashion are more greatly affected by the smallest grains. Hence, the top right panel does 
not greatly resemble the lower right panel in Figure (TUl 



4.2. Semi- Analytic Treatment 



We can further develop these ideas with a simple semi- analytic treatment. For a given 
planet mass and semi-major axis, the contrast function (see Equation H]) becomes C{s), 
where s is the particle size. We approximate the contrast function in Equation H] with the 
piecewise function 



C{s) 



for s < si 



for Si < s < S2 



(7) 



Clargo for S> S2, 

where Ciarge = 1 + Pi is the contrast for the largest particles, m is the logarithmic slope 
of the contrast in the transition regime, and Si and S2 are the particle sizes that mark the 
beginning and end of the transition regime, respectively. We fit our contrast data with this 
piecewise function and obtained the following power law estimates assuming silicate grains 
(p ~ 2 g cm~^): 



Clarge;AA,IE ~ 1 + 4.38 



0.19 
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1 Au) 



ap x-0.5 



1 AU 



(9) 
(10) 

(11) 



In the same manner as Section 13.21 we defined the contrast of any ring structure as the 
surface density within the ring, Unng, divided by the background surface density, (Tbg- The 
contrast in optical depth of a cloud containing several components of various sized particles, 
labeled with the index i, is 

(a) = -, (12) 

i 

For a coUisionless cloud with a continuous distribution of grain sizes, the contrast in optical 
depth of the composite cloud is given by 



(a) 



s^-"C(s) ds 



3^-" ds 



(13) 



where we have explicitly included the particle cross section {A{s) oc s^) and background 
surface density (o"bg('5) oc s^~"). For a coUisionless cloud, the background surface density is 
enhanced by a factor of s due to the PR time scaling as s (see Section H?T|) . and a factor of 
which describes the assumed crushing law. 

Using Equations [3, we can now integrate Equation [13] directly. Assuming Smin < Si < 
S2 < Smax and (smin/smax) < 1 whcu « 7^ 4, wc find 



large 



S2 



4— a 



_ 4-a / £2 

'-^large A-a+m Ui 



+ 



A—a+m \ s 



4-a 



iCr) = I 



In(^) + aargcln'^ 



S2 



+ m 



£2 \ _ ]^ 



1 + 



S2 
5min 



4-0 



4— a / 32 

'-"large 4-a+m [ si 



m I si 
A-a+m \ Smin 



for a < 4 
/Inf^jaai) for a = 4 

for a;>4, Q;7^4-|-m 



A— a 



1 ~l~ C*large I „ . 



A — OL 



A-a 



for a > 4, a = 4 + m. 

(14) 
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For cases in which the maximum particle size in a disk is less than S2 (see Equation [TTj) . 
simply replace all instances of S2 in Equations UM with Smax- Similarly, for cases in which the 
minimum particle size in a disk is greater than si (see Equation [TU]) . replace all instances of 

Si with Smin- 

Equations [1 
contrast in terms of Mp, Op, s 



together with Equations [HHTTl give analytic expressions for optical depth 
and Sinax- Although Equations [14] addres s all possible 
scena rios, the most plausible scenarios have crushing laws with a < 4 (e.g. iDurda et al. 
20071 ). With this assumption. Equations [H] combined with Equations [5HTT] gives 



{Ct;AA,Ie) 



1 + 4.4M^-i9 + 



(10M;~i- 



0.5 



- (1 + 4.4M*'») (l50M;-' '5a-''=) 



4-a 



for Sr 



> S2 
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Xs 



max 



where Mp 



Mp 



o' 



+ (1 - X) Smfx 

= (^^) s' -- 

V 1 AU/ ' max 



for Sr 



< S2 



1 /jm 



and X 



4-a+0.6A/^0-i8 • 



If, as in our first composite model in Section 14.11 we assurne that the particles are 
released from their parent bodies in accordance with the iDohn anvil (119691 ) crushing law and 
then spiral inward without colliding, a = 3.4. In this case. Equations UM give a contrast in 
optical depth of {Cr) ~ Ciarge in the limit Smax ^ S2- This result confirms our numerical 
results for our first composite model, shown in the upper left panel of Figure [T0| the contrast 
in optical depth is dominated by the large particles. 

For our second composite model, we forced the background density to obey a Dohnanyi 
distribution at ~ 3ap, i.e. <7bg{s) oc s~^'^, so that a = 4.5. This technique essentially 
removes the factor of s in the background surface density that results from the PR time 
scaling as s. For the composite cloud shown in Figure [101 Mp = 2.0 M^, and Op = 6.0 
AU, for which m ~ 0.68, Si ~ 1.9 /xm, S2 ~ 24 /xm, and Ciarge,AA,iE ~ 6. We let each 
simulated particle size represent a range of particle sizes using the midpoint method, which 
gives Smin ~ 0.7 /im. With these values. Equations [TU give a contrast in optical depth of 
(C't,aa,ie) ~ 2.4, in agreement with the measured contrast in the top right panel of Figure 

m 

Figure [TT] illustrates in general how a distribution of particle sizes affects the contrast of 
a ring structure. This figure compares the contrast of a collisionless multi-particle-size cloud 
(Equations [H]) to that of a single-particle-size cloud as a function of y/a^/Pmin assuming a 
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Dohnanyil (119691 ) crushing law. Both kinds of clouds have the same contrast in the adiabatic 
limit (large y/a^/ (3min), but the contribution of the smaller grains reduces the contrast else- 
where, effectively broadening the transition between the no-trapping regime and saturation 
regime. Crushing laws with a < 3.4 result in contrast curves that more closely resemble the 
single-particle-size contrast curves shown in Figure [TTl 

In a real zodiacal cloud, collisions affect the distribution of grains, even far from the 
source of the grains. Our composite dust cloud models do not include collisions and become 
unreliable for particles with coUisional times less than their PR times. Our composite models 
also lack the structural re sults of coUisi onal effects, such as the loss of particles as a function 
of circumstellar distance (IWyattI l2005l ) and any potential morphological effects in the ring 
structure. 

More sophisticated models ma y be requir e d to i nvestigate these ph enomena. However , 
since dust produced according to a iDohnanyil (119691 ) crushing law or a iDurda et al.l (120071 ) 
crushing law yields a cloud dominated by the largest grains, as we showed above, we hypoth- 
esize that resonant rings in exozodiacal clouds may often be dominated by a single particle 
size whose PR time is roughly equal to the its coUisional time. 



4.3. Ring Detectability 

The detectability of a resonant ring structure depends on many factors specific to the 
telescope being used and the observing conditions. We address this complicated issue by 
imposing one simplifying assumption: a minimum detectable optical depth ring contrast of 
1.5. This assumption likely underestimates the sensitivity of a TPF-like mission to rings in 
exozodiacal clouds analogous to the solar zodiacal cloud. In such a cloud, a ring 0.4 AU 
wide located at 1 AU from the star has ~ 15 times the total flux of an Earth-like planet at 1 
AU, even for a contrast of unity. Our assumption, conservative on the basis of photon noise 
alone, allows for the possibility of unknown systematic noise that could hinder the detection 
of extended structures. 



Figure [12] shows the minimum detectable planet mass a s a function of se mi- major axis 
and maximum dust particle size based on Equations [15] and a iDohnanyil (11969! ) crushing law. 
The masses and semi-major axes of Ea rth, Mars, and the p lanet OGLE-2005-BLG-390Lb, 
detected by the microlensing technique ( Beaulieu et al.ll2006l ). are marked for reference. This 
plot shows that an Earth-mass planet at 1 AU might be detectable if the ring contains grains 
more than a few tens of microns in size and a planet with mass equal to a few times that 
of Mars might be detectable near 10 AU if the ring contains grains more than one hundred 
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microns m size. 

The detectability of a ring structure depends upon the size distribution of d ust within 



the rin g structure. Dust produced according to a crushing law less steep than the iDohnanyi 



(119691 ) crushing law {a < 3.4) will result in more highly contrasted ring structures because 
of the increased relative contribution of the large grains. For crushing laws with a = 3 and 
a = 2, the curves of constant maximum particle size shown in Figure [12] shift downward by 
a factor of approximately 1.25 and 1.55, respectively. 

These values are subject to the assumptions of our simulations, which do not include a 
dynamically hot component in the dust cloud. This component would reduce the contrast 
in the ring, making planets harder to detect for a given cloud mass. So the detection limits 
shown in Figure [T^ should be thought of as best-case scenarios. 



5. Caveats 

Our simulations include a number of simplifying assumptions, which we summarize here. 
We ignored the effects of dynamically hot dust, like dust that might come from comets. 
Trapping probability decreases dramatically for particles on highly eccentric and inclined 
orbits, so we expect dynamically cold dust to dominate any resonant debris disk structure. 
As a first approximation, we can treat the contribution from the dynamically hot dust as 
a constant surface density cloud component, which reduces the contrast of any structure 
formed from the dynamically cold component. Estimates of the ratio of astero idal dust to 



cometary dust in our solar system range from 1:10 to 7:10 (llpatov et al.l 120071 ). For other 
systems, this ratio is also unknown. 

Our simulations also assumed a single planet on a circular orbit around a Sun-like star. 
We have performed trial simulations of our solar system and demonstrated that the presence 
of Jupiter may reduce the Earth's ring contrast. Other multiple-planet systems may also 
exhibit a similar effect. Additionally, planets on ecce ntric orbits give rise to add itional MMRs 



with different capture probabilities and geometries (IKuchner fc HolmanI 120031 ). 



The ring contrasts of inclined systems can vary significantly depending on the inclination 
and radial extent of the dust cloud. In edge-on systems, resonant features can overlap as 
seen from the Earth, complicating their interpretation. The contrasts we provide are useful 
only to systems for which projection effects can be taken into account. 



Finally, our multi-particle-size models demonstrate the subtlety of collisional effects in 
dust clouds. Collisional effects can determine the relative populations of large and small 
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grains and potentially alter the morphology of the ring structures. Our simulations can not 
yet handle these effects in detail. 



6. Conclusions 

We have implemented our own hybrid symplectic integrator for the ri-body problem and 
used it to simulate coUisionless debris disks, taking into account solar wind and drag effects. 
Each simulation contained 5,000 particles. We found that this number of particles suffices 
to populate the dominant MMRs of a low-mass planet with an accuracy at the few percent 
level, yielding for the first time models of the surface brightness distributions of exozodiacal 
clouds that we can use to quantitatively study the contrasts of resonant features — not just 
their geometries. 

We generated a catalog of resonant structures induced by a single planet on a circular or- 



bit around a Sun-like star, available online at http:/ /asd. gsfc.nasa.gov/Christopher.Stark/catalog.php 



We investigated 120 sets of model parameters, spanning a range of planet masses, planet 
semi-major axes, and values for /3, assuming dust grains launched from orbits with low Cdust 
and Zdust- The resulting ring structures exhibited leading-trailing asymmetries, gaps near the 
locations of the planets, and sharp inner edges at ~ 0.83ap. 

We performed a detailed analysis of the surface density contrasts of the rings (Figure 
|9]). We showed that for a planet on a circular orbit, the contrast and morphology of the 
rings are to good approximation functions of only two parameters, Mp and y/a^//3, for a 
given stellar mass and distribution of dust sources for Mp < 5 and j3 < 0.25. Equation H] 
summarizes the contrasts of our single-particle-size models as a function of these parameters. 
Considering only the dynamically cold particles analogous to particles released by asteroids 
in the solar system, we find that terrestrial-mass planets are capable of producing resonant 
ring structures with azimuthally averaged contrasts up to ~ 7 : 1. 

By combining our simulations of grains with particular f3 values, we assemb led multi 



particle-size models of 25,000 particles each. Releasing the particles according to a lDohnanyi 



(1l969l ) crushing law without any subsequent collisional processing results in composite clouds 
whose optical depths are dominated by large particles; large particles will dominate images 
of these clouds in visible light and throughout the IR. Based on these composite models, 
we suggested that the best current models for exozodiacal clouds are those with a narrow 
range of grain sizes corresponding to grains whose collision time roughly equals their PR 
time. Future models should account for processes like grain-grain collisions that destroy 
large grains. 



- 19 - 



Equations [M] and [15] provide semi-analytic predictions for the contrast in optical depth 
of a multi-particle-size cloud of dynamically cold grains. For ring structures composed of 
silicate grains released according to a Dohnanyi crushing law {a = 3.4), Equation [T^ gives 
an approximate contrast of 



1 _L 4 4 ( JHe. 



0.19 



-52 



1 fim 



-0.6 



-0.3 / M, 



M, 



-0.57 



for Smax > S2 
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-0.3 



for Sr 
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(16) 

where (CrjAA.m) is the ratio of the azimuthally averaged optical depth in the ring structure 
to the azimuthally averaged background optical depth, Smax is the maximum grain size in 
the ring structure, and S2 is given by Equation [TTJ For the case Smax > -52, the first two 
terms in Equation [T6l represent the contrast in the adiabatic limit. The remaining term (and 
the terms in the Smax < S2 case) represents deviations from this limit for smaller particles or 
smaller semi-major axes. 

We plotted the mass of the smallest planet that could be detected through observation 
of a resonant ring structure as a function of planet semi-major axis and parti cle size in 
Figure [T2l We assumed a cloud composed of a range of particle sizes adhering to a iDohnanyi 
fll969[ ) crushing law and a minimum detectable optical depth contrast of 1.5:1. We found 
that planets with masses just a fraction of the Earth's may form detectable ring structures 
if the rings harbor grains more than several tens of microns in size. 
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Fig. 1. — Maximum fractional error in energy during a 3,000-year integration as a function of 
perihelion distance for two scenarios: a two-body system of the Sun & Jupiter (solid line) and 
a three-body system of the Sun, Jupiter & Saturn (dashed hne). For the two-body system, 
Jupiter's perihelion distance is plotted. For the three-body system, Saturn's perihehon 
distance was altered while Jupit er's rern ained fixed. The inclinations and eccentricities of 
both planets remained fixed, cf. IDLL98I Fig. 3. 
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Fig. 2. — Top: eccentricity as a function of time for a dust particle with (3 = 0.2 and 
s w = 0.3 5. Bottom: s emi-major axis as a function of time. Our resu lts match those of 



MMMO2I (cf. IMMMOI Fig. 1) and agree with the analytical solution (IWvatt fc Whipple 
1950h. 
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Fig. 3. — The integrated trajectory of comet P/Oterma during a close encounter with Jupiter 
as viewed in the frame centered on and rotating with Jupiter with the Sun on the negative 
X-axis. Shown are the results of a Buhrsch-Stoer integrator and our hybrid symplectic 
integrator for four values of integration time step. The hybrid symplectic resufrs overlap the 
Buhrsch-Stoer resufrs for a timestep of 1 day (cf. IC99l . Fig. 4). 
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Fig. 4. — Comparison of our hybrid symplectic integrator with a Bulirsch-Stocr integrator. 
Left: surface density histogram for 1,000 particles in the Sun-Earth system using a Bulirsch- 
Stoer integrator. Middle: surface density histrogram for the same initial conditions using our 
hybrid symplectic integrator. Right: Bulirsch-Stoer histogram minus the hybrid symplectic 
histogram (image is in units of a, the \/n Poisson noise associated with the histograms). 
Except in a handful of pixels, the difference is roughly consistent with Poisson noise. 
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Fig. 5. — Population of Neptune's MMRs for three independent simulations of 5,000 particles 
each (shown in green, red, and black). Populations of the 2:1 and 3:2 resonances differ among 



the three simulations by 6.4% and 4.3%, respectively (c.f. iMoro-Martfn fc Malhotral 12002 
Fig. 5). 
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Fig. 6. — Surface density distributions for four of the 120 simulations (scale is relative). 
The star is located at the center of the image and the planet is marked with a white dot. 
The planet orbits counter-clockwise in these images. Integrations were truncated at half the 
planet's semi-major axis. The simulations shown on the right have different values of Op and 
P, but the same value of ya^//3. Their surface density distributions are nearly identical; 
their difference is consistent with Poisson noise (see Section [32]) ■ 
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Fig. 7. — Left: the angular size of the "gap" in the ring structure around the location of the 
planet in our simulations versus the contrast of the ring structure. Right: the sine of the 

1 /2 

prograde shift plotted against the function — (3)/{Mpap ). We removed all data with 
C'aaje < 1-6 from these plots. Solid hnes show hnear fits to the data. 
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Fig. 8. — The azimuthally-averaged contrast measured at the inner edge of the ring structure 
(see Section 13.21 for definition of contrast) as a function of (3 (left panel) and planet mass 
(right panel). Both figures show a transition from a no-trapping regime to a saturation 
regime where contrast is independent of semi-major axis. 
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Fig. 9. — The contrast in surface density of the ring structure compared to the background 
cloud for all combinations of Mp, ap, and (3 (see Section [3^ for definitions of contrast). The 
contrast is only a function of two parameters: planet mass and y^/p. The solid lines are 
fits to the data (see Equation Hj) . 
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Fig. 10. — Comparison of the optical depths for a composite cloud formed by two different 
methods for Mp = 2.0 and Op = 6.0 AU. The planet, marked with a white dot, orbits 
counter-clockwise in these images. Top-left: A composite coUisionless clou d where the part i- 
cles are released with a size distribution equal to the crushing law used by iDohnanyil (119691 ). 
Top-right: The same composite cloud, but formed by forcing the surface density outside of 
the ring structure to obey a Dohnanyi distribution. Bottom-left: The optical depth of the 
smallest particles included in the composite clouds. Bottom-right: The optical depth of the 
largest particles included in the composite clouds. The largest particles dominate the optical 
depth in a cloud of particles released with a Dohnanyi crushing law. 
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Fig. 11. — Contrast in optical depth for multi-parti cle-size clouds (s olid lines) compared to 
single-particle-size clouds (dashed lines) assuming a iDohnanyil (119691 ) crushing law (a = 3.4; 
see Equations [T6|l . From top to bottom, the six solid lines and six dashed lines correspond 
to six values of planet mass: 5, 2, 1, 0.5, 0.25, and 0.1 Mq. The contributions of the small 
grains reduce the contrasts of the multi-particle-size clouds compared to single-particle-size 
clouds with the same minimum value of (3. 
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Fig. 12. — Minimum detectable planet mass in a multi-particle-size coUisionless cloud as a 
function of semi-major axis and maximum grain size, assuming a Sun-like star, a minimum 
detect able ring contrast of Cr,AA,iE = 1-5, and dust produced according to a iDohnanyi 
( 1l969l ) crushing law {a = 3.4). Earth-like and Mars-like planets are denoted with an E 
and M, respectively. The 5.5 M® exoplanet OGLE-2005-BLG-390Lb is denoted with an O 
( jBeaulieu et al.ll2006l ). Listed values for maximum dust size in the ring structure assume 
perfectly absorbing spherical grains with mean density p = 2.0 gm cm~^ and radius Smax- 
The bold line shows the case of the solar zodiacal cloud, for which the observed emission is 



dominated by 30 fim grains (iFixsen fc Dwekll2002l ). The dashed lines show typical inner and 
outer detection limits for a mission similar to TPF. 



